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Abstract 

In this paper, we propose a new technique for two-dimensional phase unwrapping. The unwrapped 
phase is found as the solution of an inverse problem that consists in the minimization of an energy 
functional. The latter includes a weighted data-fidelity term that favors sparsity in the error between the 
true and wrapped phase differences, as well as a regularizer based on higher-order total-variation. One 
desirable feature of our method is its rotation invariance, which allows it to unwrap a much larger class 
of images compared to the state of the art. We demonstrate the effectiveness of our method through 
several experiments on simulated and real data obtained through the tomographic phase microscope. 
The proposed method can enhance the applicability and outreach of techniques that rely on quantitative 
phase evaluation. 
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1 Introduction 


Two-dimensional phase unwrapping is an essential component in a majority of techniques used for quantita¬ 


tive phase imaging. A typical example of standard application is tomographic phase microscopy CFYB+07 


where the phase of the transmitted wave-field must be first unwrapped before being interpreted as a line 
integral of the refractive index along the direction of propagatio n. Phase unwrapping is also used for estima¬ 
tion of terrain elevation in synthetic aperture radar RHJ+00 , wave/fat separation in magnetic resonance 


imaging |SNPG95 , and estimation of wave-front distortion in adaptive optics Fri77 . Accordingly, improve¬ 


ments in phase unwrapping methodology can enhance the applicability and outreach of techniques that rely 
on quantitative phase evaluation. 

The centrality of phase unwrapping has resulted in the development of many practical solutions to the 
problem. Such solutions include direct approaches based on path-following such as the two-dimensional exten¬ 
sion of the well known Itoh’s method Ito82 as well as more evolved strategies based on branch cuts GZW88 


or quality maps HBLG02 . The current trend in the literature is to formulate the task of phase unwrapping 


as an inverse problem in a path-independent way. Earlier works have proposed to minimize the quadratic 
between the true and wrapped phase differences [Hun79 TT88 GR94 . Marroquin and Rivera |MR95 


error 


have applied Tikhonov regularization to improve and stabilize the performance of the least-squares approach. 
They showed that the intro duction o f the regularization term permits the algorithm to cope with noise and 
missing data. Huang et al HTZ+12 showed that the performance of phase unwrapping is further improved 


when replacing Tikhonov with total-variation regularization. In the context of phase unwrapping, Ghiglia 
and Romero GR96 recognized the tendency of quadratic-penalty to smooth image edges and proposed to 


minimize a more general ip norm-based criterion as an alternative. They found that the performance of 
phase unwrapping improves when 0 < p < 1, albeit at the increase in computational cost. Similarly, Rivera 
and Marroquin RM04 have investigated nonconvex optimization strategies relying on half-quadratic reg¬ 
ularization. More recently, Gonzalez and Jacques GJ14 have proposed an iterative unwrapping method 


based on mi nimization that additionally promotes sparsity of the solution in the wavelet domain. Ying 
et al. YLJ+06 have proposed an iterative method based on dynamic programming that models the phase 
as a Markov random field. Bioucas-Dias and Valadao BDVQ7| have proposed an energy functional based 
on generalized ip norm and corresponding minimization algorithm that relies on graph-cut methods. Their 
PUMA algorithm in its original form and its noise tolerant extensions 
considered state of the art. Mei et al 


BDKAE08 VBD09 are currently 


MKCG13 have proposed an application specific method that jointly 


unwraps and denoises time-of-flight phase images using a message-passing algorithm. An extended review 
of this topic, along with related algorithmic ideas, can be found in the book 


GP98 and tutorial Yin06 


In this paper, we propose a new variational-reconstruction approach for phase unwrapping that is robust 
to noise. In particular, our aim is to improve on state of the art by introducing an improved energy functional 
and demonstrating its benefits. The main contributions of this paper can be summarized as follows: 

• Eormulation of the unwrapping as an optimization problem where the data-fidelity term penalizes the 
weighted ^i-norm of the error in a way that is invariant to rotations. Our formulation thus allows 
the phase image to contain edges that are of arbitrary orientation, which is distinct from traditional 
approaches in literature GR96[ RMQ4, YLJ+Q6 BDV07 GJ14 . 


Use of a non-quadratic regularization term that allows our method to cope with noise, while still 
preserving sharp edges in the phase image. Our regularizer consists of a higher-order extension of 
total-variation (TV) that is currently considered state of the art in the context of resolution of linear 
problems in biomedical imaging [LWU13 . 


inverse 


Design of a novel iterative algorithm for phase unwrapping. The algorithm approximates the minimum 
^o-norm solution by solving a sequence of weighted -^i-norm minimization problems, where the weights 
at the next iteration are computed from the value of the current solution. Since our energy functional 


is non-smooth, we rely on a well known alternating direction method of multipliers (ADMM) BPC+11 
to decompose the minimization into a sequence of simpler operations. 
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This paper is organized as follows. In Section we introduce our formulation of the phase unwrap¬ 
ping, and discuss the relevance of this new approach for obtaining high-quality solutions in practice. In 
Section we derive our reconstruction algorithm. In Section we conduct experiments on simulated and 
real phase unwrapping problems, and compare our method with the state of the art from both qualitative 
and quantitative standpoints. We summarize and conclude our work in Section 


2 Problem formulation 

We consider the following observation model 


= W (0) = 4 > — 27r/c 


( 1 ) 


where k G and ^ The vectors -0 and 0 represent vectorized versions of the wrapped and 

unwrapped phase images, respectively. The wrapping is represented by a component-wise function W that 
is defined as 


>V(0) = [((0-h7r) mod 27r) — tt] 


[-7r,7r). 


When noise is part of the measurements, we assume that the unwrapped phase vector (j) in our model 
represents the noisy version of the true phase Lp. 

Generally, the two-dimensional phase unwrapping problem is ill-posed. However, it can be solved exactly 
in the noi seless scenario, when the phase 0 satisfies the two-dimensional extension of Itoh’s continuity 
condition |lto82 . Let D : ^ denote the discrete counterpart of the gradient operator and let 


D0 = 


^x4> 

Dy(/) 


( 2 ) 


where and denote the finite-difference operator along the horizontal and vertical directions, respec¬ 
tively. If, for a given pixel n G [1,..., A^], the unwrapped phase 4> satisfies 


= ^([D,(/>]„)2 + ([D,0]„)2<7r, 


then, we have the equality 


[D(/>]„ = >V([DV’]„). 


(3) 


(4) 


Here, \Dy4>]n) denotes the n-th component of the gradient D0. Relation 0 suggests that 

two-dimensional phase unwrapping may be accomplished by a simple phase summation, provided that 
is satisfied at all pixels n. Note that the formulation in ^ imposes the Itoh’s continuity condition on both 
gradient components simultaneously due to the norm inequality 


max (|a| , |6|) < \/a^ + 6^, 


which holds for any a, 6 G M. 

In practice, however, condition will not be fulfilled at all pixel locations due to the presence of sharp 
edges and of measurement noise. Yet, it can still be expected to hold for the great majority of pixels of the 
unwrapped phase image. We thus formulate phase unwrapping as the following minimization problem 

0 = arg min {T>(0) + r7^(0)} , (5) 

where V is the data-fidelity term and IZ is the regularization term, to be discussed shortly. The convex 
set ^ = {0 G C enforces the first pixel of the solution 0 to match the first pixel of 

the wrapped phase -0, which removes the additive constant ambiguity present in phase unwrapping. The 
parameter r > 0 controls the amount of regularization. 
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The data-fidelity term in (§ is given by 


N 




(6) 


n=l 


where w G are positive weights. It is intended to relax the strict equality 

In the unweighted case, i.e. when Wn = 1 for all n G [1,..., A/"], V corresponds to ^i-norm penalty on 
the magnitudes of e = Y>(f) — >V(D'0). It can be interpreted as a convex relaxation of -^o-norm penalty that 
enforces sparse magnitudes of e. This implies that our data-term favors (j) whose gradient T>(j) agrees with 
>V(D'0) on most of the pixels. Moreover, our data-fidelity term V penalizes both horizontal and vertical 
components of e in a joint fashion. This is significantly different from traditional formulations in the 


literature GR96[[GJ14 , where the ^i-norm is penalized in a separable fashion as ||€||^^ = 
fact, there is a clear analogy between our formulation and the isotropic^ i.e. rotation invariant, form of 


total-variation (TV) that is often used for edge-preserving image restoration BT09 . Similarly, the separable 
-^i-penalty ||e||^^ is analogous to the anisotropic form of TV. The arbitrary orientation of edges in a typical 
image makes isotropic TV penalty a preferred choice for image restoration. The numerical experiments 
presented in Sec. illustrate that indeed our method based on isotropic formulation can unwrap a larger 
class of images compared to other state-of-the-art phase unwrapping methods. 

It has been reported in several works GR96 BDV07 that nonconvex approaches based on ^^-norm 


penalization of e, with 0 < p < 1, further improve the performance of phase unwrapping. In particular, 
-^o-norm is generally accepted as the most desirable in practice. One of the properties of -^i-norm that 
distinguishes it from Iq is that it takes into account the actual values of the magnitudes of e, whereas Iq- 
norm disregards this information and only counts the support. One possible strategy for selecting weights in 
the context of £i minimization proposed by Candes et al. GWB08 is to pick w such that it counteracts the 


influence of the magnitude on the -^i-norm. For example, suppose that the weights were inversely proportional 
to the true magnitudes 

if ll[e]nll ^2 = 0. 


= 


(7) 


Then, if there are exactly m pixels violating Q, the minimizer in (§ is guaranteed to find the solution 
corresponding to £o data-fidelity term. The large values in w force the solution 0 to concentrate on the 
pixels where the weights are small, and by construction these correspond precisely to the pixels where the 
magnitudes of e are nonzero. Typically, these pixels correspond to the area of the phase image that contains 
a sharp edge. It is clearly impossible to construct the precise weights 0 without knowing the unwrapped 
phase (j) itself, but this suggests more generally that large weights could be used to discourage nonzero 
magnitudes in e, while small weights could be used to encourage nonzero magnitudes in e. 

As regularization term in (|^, we propose to use the Schatten norms of the Hessian matrix at every pixel 
of the image LWU13 . Specifically, we set 


7^(0) 4 ||H0||i,< 


N 


4 ^(cTi([H</)]„) + a2([H</)]„)), 


(8a) 

(8b) 


n=l 


where H : ^ ^Nx2x2 discrete Hessian operator and is the nuclear norm that is computed 

by summing singular values ai{[tlcj)]n) and <72([H(/)]^) of the Hessian matrix at position n. There are four 
major advantages of using such Hessian Schatten-norm (HS) regularization: 


• As a higher order regularizer HS avoids the staircase effect of TV and results in piecewise-smooth 
variations of intensity in the reconstructed phase image. Accordingly, this makes HS particularly well 
suited for biological and medical specimens that consist of complicated structures such as filaments. 
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Figure 1: Phase unwrapping of 128 x 128 image of a truncated Gaussian function, (a) true phase, (b) 
wrapped phase, (c) Goldstein’s algorithm (GA), (d) least-squares (LS), (e) iteratively reweighted LS (IRLS), 
(f) PUMA, (g) proposed method with uniform weights, (h) proposed method with adaptive weights (IRTV). 


• Similarly to the data-fidelity term P, our regularizer IZ is convex and rotation invariant |LWU13 . 
Rotation invariance implies that we can expect it to work equally well on phase images that contain 
objects of arbitrary orientations. 


• HS regularization has been shown to be state of the art in resolution of linear inverse problems. In 

that HS consistently outperforms other popular regularizers such 


particular, it was shown in 
as Tikhonov, wavelet, and TV. 


LWU13 


• Gonvexity and its algebraic structure make HS amenable to efficient algorithmic implementation and 
thus practical for large scale inverse problems that are typical in imaging. 


To demonstrate the performance of our variational approach ([^, we present a phase unwrapping exper¬ 
iment on a synthetic image consisting of a 2D Gaussian function of amplitude 12 and standard deviation 
20 that has been truncated along a line of arbitrary orientation (here about 40^). Fig. 
results for four standard unwrapping methods such as Goldstein’s algorithm (GA) [GZW88] 

iteratively reweighted LS (IRLS) with data-dependent weights that approximate the ^o-norm 
BDVQ7| . We additionally illustrate the performance of our method with the 


(LS) [GRM 


I: 


illustrates the 
least-squares 


penalty GR96 , and PUMA 


emphasis on the influence of weights w over the final solution 0. Accordingly, we show the solution of ^ 
with uniform and data-dependent weights. As expected, all algorithms perform equally well in the contin¬ 
uous region of the image. On the other hand, our approach is the only one that accurately captures the 
discontinuous region of the unwrapped image. This is expected due to rotation invariance of our energy 
functional. Additionally, we note that the LS method, which is also based on rotation invariant energy 
functional, fails to preserve the edge due to excessive smoothing. Finally, a careful inspection of Figs. 0(g) 
and (h) reveals that the edge is much sharper when the weights w are selected in a data-dependent fashion 
as explained next in Sec. 


3 Reconstruction algorithm 

We now describe our computational approach based on the convex optimization problem (|^. The iterative 
scheme alternates between estimating 0 and redefining the weights w as follows. 

1. Initialization: Set iteration number to t = 1. Select an initial phase (pP G and set = 1 for each 

n = 1,..., V. 
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Figure 2: Set of standard test images. From left to right: Cameraman, Peppers, Lena, Man, Lake, and 
Barbara. 


2. Optimization: For a fixed compute the phase image cj)^ by solving (§. Also, compute the 

auxiliary variable — >V(D'0). 

3. Weight adaptation: For each n = 1,..., A" 

if ^min ^ II 11^2 — ^max 

if ||[e*]„|k, >e„,ax (9) 

if ||[e*]„||£2 < eniin 


W„ = < 


ll[e‘]nlU2 


tmax 

1 


4. Stop on convergence or when t attains a specific maximum number of iterations tmax- Otherwise, 
increment t and proceed to step 2. 

The two parameters emin and emax in step 3 provide stability and avoid divisions by zero. For our experiments 
in Sec[^ we set e^ax = 1/emin = 10. The optimization in step 2 will be discussed shortly. 

Although the initial phase can be set to an arbitrary vector in in practice, a warm initialization 
leads to a smaller number of iterations required for convergence, and hence faster unwrapping times. In our 
experiments, we found that the solution of LS can serve as a computationally inexpensive way of obtaining 
a good initialization. 

Using an adaptive approach to construct the weights progressively improves the unwrapping around the 
discontinuities in the phase image. These phase discontinuities might be due to a presence of a sharp edge 
or due to strong noise. Even though the early phase estimate may be inaccurate, the largest coefficients of e 
are most likely to be identified with a phase discontinuity. Once these locations are identified, their influence 
is downweighted in order to gain in sensitivity for identifying the remaining regions of the phase image. 

The step 2 of the algorithm requires the resolution of the non-smooth optimization problem (§. We 
perform this minimization by designing an augmented-Lagrangian (AL) scheme NW06 . Specifically, we 
seek the critical points of the following AL 


N 


£((/),€,s) = +r7^((/)) 


n=l 


+ (e _ D</) + d) + pe - D</. + 


€ — D0 d — 
P 


t‘2 


N 




n=l 


where d = yV(D'ip) is the data vector, s G is the dual variable that imposes the constraint e = D0 —d, 

and p > 0 is the quadratic penalty parameter. Traditionally, an AL scheme solves the problem (§ by 
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Table 1: Performance of the five considered phase unwrapping methods in terms of SNR on the set of 
standard test images for various amplitudes 


Img/Amp 

SNR (dB) 

Img/Amp 

SNR (dB) 

GA 

LS 

IRLS 

PUMA 

IRTV 

GA 

LS 

IRLS 

PUMA 

IRTV 

Cameraman 

4 

33.99 

38.76 

32.74 

35.72 

oo 

Man 

4 

oo 

oo 

oo 

oo 

oo 

5 

23.38 

25.79 

24.07 

23.89 

25.65 

5 

26.29 

32.31 

oo 

oo 

oo 

6 

14.10 

15.51 

19.06 

17.41 

19.98 

6 

10.12 

24.55 

30.47 

29.43 

oo 

7 

-1.31 

7.57 

15.64 

12.73 

16.09 

7 

11.70 

18.14 

24.23 

27.46 

29.22 

8 

0.00 

2.35 

0.38 

0.73 

0.92 

8 

8.75 

10.17 

21.89 

12.08 

25.26 

9 

1.63 

2.64 

2.02 

1.90 

2.17 

9 

6.95 

6.75 

12.69 

8.36 

13.23 

Peppers 

4 

oo 

oo 

oo 

oo 

oo 

Lake 

4 

oo 

oo 

oo 

oo 

oo 

5 

36.69 

oo 

oo 

oo 

oo 

5 

oo 

oo 

oo 

oo 

oo 

6 

23.64 

22.60 

oo 

oo 

oo 

6 

40.10 

oo 

oo 

oo 

oo 

7 

14.40 

11.67 

oo 

26.39 

oo 

7 

26.62 

25.21 

oo 

41.44 

oo 

8 

7.27 

8.23 

23.96 

14.93 

27.62 

8 

15.99 

13.84 

18.02 

17.69 

18.10 

9 

5.58 

3.97 

15.69 

15.64 

15.75 

9 

-0.88 

2.84 

15.91 

15.33 

16.52 

Lena 

4 

oo 

oo 

oo 

oo 

oo 

Barbara 

4 

oo 

oo 

oo 

oo 

oo 

5 

oo 

oo 

oo 

oo 

oo 

5 

oo 

oo 

oo 

oo 

oo 

6 

38.56 

27.77 

oo 

oo 

oo 

6 

oo 

oo 

oo 

oo 

oo 

7 

23.62 

18.70 

oo 

oo 

oo 

7 

43.54 

40.53 

oo 

oo 

oo 

8 

12.20 

13.33 

36.29 

30.30 

oo 

8 

35.60 

30.07 

44.70 

oo 

oo 

9 

-2.97 

7.33 

25.50 

15.51 

30.47 

9 

31.01 

22.01 

oo 

45.72 

oo 


Table 2: Phase unwrapping results for the noisy Man image of size 128 x 128 and amplitude a = 6. 


Method / Input SNR 

16 dB 

18 dB 

20 dB 

GA 

14.40 

19.03 

18.72 

LS 

23.57 

25.75 

27.97 

PUMA 

31.13 

31.96 

34.96 

IRTV (t = 0) 

28.13 

30.21 

34.96 

IRTV (t = 10-3) 

28.52 

31.97 

oo 

IRTV (t = 10-2) 

35.02 

34.98 

oo 

IRTV (t = 10-1) 

18.86 

24.56 

28.46 


alternating between a joint minimization step and an update step as 


(0^,€^)^ argmin {£(0,€,s^ ^)} 

(lOa) 

sk ^ 5^-1 4-p(e^ + d). 

(10b) 


However, 

separate 

direction 


the joint minimization 
(10a) into a succession 
method of multipliers 


step can be computationally intensive. To circumvent this problem, we 
of simpler steps. This form of separation is commonly known as alternating 
(ADMM) BPC+ll| and can be described as follows 


0^ ^ argmin {£(0, ^)} (11a) 

^ argmin{£(0^€,s^-^)} (11b) 
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Figure 3: Phase unwrapping results for Lena image of size 256 x 256 when the amplitude is a = 9. (a) LS 
(SNR = 7.33 dB); (b) IRLS (SNR = 25.50 dB); (c) PUMA (SNR = 15.51 dB); (d) IRTV (SNR = 30.47 dB). 


By ignoring the terms that do not depend on 0, the step (11a) can be expressed as 


cj)^ ^ arg min 




- 7 ^( 0 ) 

P 


( 12 ) 


with + d + s^~^/p. This step corresponds to a classical Hessian Schatten-regularized linear 

inverse problem. The solution of this problem can be efficiently solved with the publicly available software 
that has been described in [LWU13| . Similarly, the step (11b) can be simplified as follows 


k ' \ P w k 

e ^ arg mm < - \\^ — y 


N 


+ '^n II [^]n 11^2 


n=l 


with = D(/)^ — d — ^ /p. This step is solved directly by component-wise application of the following 

shrinkage function 

r( 2 /;r) = argmin j^||x-i/|||^ +r||x||^,| 

= max( 112 / 11^2 -UO) n^i^. 
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Figure 4: Phase unwrapping results for Man image of size 256 x 256 when the amplitude is a = 7. (a) LS 
(SNR = 18.14 dB); (b) IRLS (SNR = 24.23 dB); (c) PUMA (SNR = 27.46 dB); (d) IRTV (SNR = 29.22 
dB). 


Thus, we can express (11b) as 


[€'']„ ^ r ( -d- s’^-^/p] ^ ; Wn/p) , 


for every n = ... ,N. 

While the theoretical convergence of our algorithm requires the full convergence of ADMM inner it¬ 
erations (11), in practice, we found that, by using a sufficiently high number of iterations /Cmax with an 


additional stopping criterion, our algorithm achieves excellent re sults as ill ustrated in Sec. In particular, 
we implemented the standard criterion suggested by Boyd et al. [BPC+ll , where ADMM is stopped when 
then primal and dual residuals are small 


||e'= - - d\U, < 

\\pB{<f>>^ - </)■ 


k-1 


^ dinner? 


(13a) 

(13b) 


where the constant dinner > 0 controls the desired inner tolerance level. In all our experiments, we set 
dinner 10 aud A^max 2000. 
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Figure 5: Phase unwrapping results for Beads image of size 512 x 512 obtained from the phase of the 
transmitted field, (a) wrapped phase; (b) GA; (c) LS; (d) IRLS; (e) PUMA; (f) IRTV. Scale bar, 5 /am. 


For outer iterations, our algorithm relies on a separate number of tmax iterations with a distinct stopping 
criterion based on measuring the relative change of the solution in two successive iterations as 


!!</>*Ik 


<5, 


— ^outer5 


(14) 


where (Pouter > 0 controls the desired outer tolerance level. In the experiments, these constants were set to 

^max — 10 and (^outer — 10 

It is important to note that the solution of our iterative method is not consistent in the sense that 
the rewrapped phase W(0) is not necessarily equal to the measured phase -0. The possible inconsistency 
of our solution comes from the fact that we are using a continuous optimization for solving an inherently 
discrete optimization problem (i.e., addition and subtraction of integer multiples of 27r). Accordingly, any 
method relying on continuous optimization such as LS, IRLS, or the method proposed here, may result 
in an inconsistent solution. Path-following methods such as Goldstein’s algorithm or discrete optimization 
algorithms such as PUMA always return consistent solutions. Gonsistency, however, is easily achieved with 
a single post-processing step that was proposed by Pritt [Pri as 

0'^ 0 + >V('0 - 0), (15) 


where (j)' is consistent with the wrapped phase -0. In the sequel, we perform a single application of the 
operator (15) to the outputs returned by the continuous optimization algorithms to make their solutions 
consistent with the measurements. 

To conclude, we described a method for minimizing the proposed objective functional. While this forces 
us to revert to a more costly iterative scheme (instead of finding the solution directly), it allows us to obtain a 
variational formulation that incorporates the most efficient ideas that have appeared for phrase unwrapping 
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and the stabilization of ill-posed inverse problems. The optimization itself is performed iteratively by relying 
on ADMM to reduce the problem to a succession of straightforward operations. The final computational 
time required to unwrap a given image depends on the total number of iterations, which in turn depends on 
the severity of wrapping. For example, it took us about 3 minutes to unwrap the Gaussian signal in Fig. 
with a MATLAB implementation of our algorithm running on an iMac using a 4 GHz Intel Gore i7 processor. 
In Sectionwe illustrate that the improvement in reconstruction quality can be rather substantial. We thus 
believe that the method should be of interest to practitioners that rely on quantitative phase evaluation. 


4 Experiments 


Based on the above developments, we report the results of our phase unwrapping method in simulated and 
practical configurations. In particular, we compare the results of our approach, which we shall denote IRTV, 
against those obtained by using four alternative methods; namely, Golstein’s algorithm (GA) [GZW88 , least- 
squares (LS) [GR94 , iterative reweighed LS approach (IRLS) with data-dependent weights that approximate 
the ^o-norm penalty [GR96 , and PUMA [BDVQ7| . An implementation of PUMA is freely available at 
http://www.lx.it.pt/~bioucas/code.htm, Mirroring examples provided by the authors of PUMA, we 
use it with a nonconvex quantized potential of exponent p = 0.5 where quadratic region threshold is set to 
0.5. As suggested in the PUMA code, we use cliques of higher order by considering 4 displacement vectors 
(1,0), (0.1), (1,1) and (-1,1). 

In the sequel, a first set of experiments with simulated phase wrapping evaluates the algorithms quan¬ 
titatively by using the true phase (f) for comparison. The experiments that follow allow to determine the 
relevance of our approach on real data in the context of tomographic phase microscopy. 


4.1 Synthetic Data 

In this set of experiments, we use a set of 6 grayscale 256 x 256 shown in Fig. [^from the USC-SIPI database 
at http: //sipi. use. edu/database/, After normalization of its amplitude to [0, a] , with a G [4, 5, 6, 7, 8 , 9], 
each image is used to generate a distinct wrapped phase image according to our observation model Q. As 
image amplitude a increases, the edges in the true phase (j) become sharper, thus making the unwrapping 
task increasingly more difficult. 

Given the wrapped phase image -0, our goal is to determine how accurately the true phase 0 can 
be reconstructed by the standard and proposed methods. As reconstruction parameters, our algorithm 
uses r = 10“^ and e^ax = 1/emin = 10 in all synthetic experiments. As mentioned before, we set the 
number of outer and inner iterations to tmax = 10 and /cmax = 2000, respectively, with additional stopping 
criteria (13) and jp). We also use 10 iterations for finding the solution of Hessian Schatten-regularized 
inverse problem (12). In principle, the positive scalar p of ADMM can either be fixed to some predeterm ined 
value or adapted according to the distance to the constraint via the scheme described in [BPG^ll . In 
all our experiments, we used adaptive p; however, we found that in practice fixed values of p = 0.1 and 
p = 1 work equally well. Iterative algorithms IRLS and IRTV were initialized with the LS solution. Upon 
convergence, the solutions of all the variational algorithms were made consistent with a single application of 
the operator ([T^. 

In Table 1 we report the signal-to-noise ratio SNR (dB) = lOlog^^g unwrapped 

phase for alFthe methods considered. When the SNR is more that 300 dB, we consider the unwrapping to 
be exact and set the corresponding value in the table to infinity. As can be seen from the table, our 
method, which is labeled IRTV, consistently provides better reconstruction results for nearly all images and 
all amplitude values. Specifically, our method successfully unwraps of Cameraman^ Lena, and Man phase 
images for which all other standard methods fail. Note that for Cameraman at amplitudes 5, 8, and 9 LS 
yields unwrapped phase images of higher quality. While for a = 5 the difference between LS and IRTV is 
modest (less that 0.15 dB), for a = 8 and a = 9 all the methods completely fail at unwrapping (SNRs below 
3 dB). 
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Beyond the SNR comparisons, the effectiveness of the proposed method can also be visually appreciated 
by inspecting the more difficult scenarios of Lena and Man images presented in Fig. and Fig. From 
these examples we can verify our initial claim: our phase unwrapping method recovers oriented sharp edges 
more accurately than other algorithms thus providing an overall boost in performance. 

In order to illustrate the robustness of IRTV to noise, we now consider a simple scenario where the 
unwrapped phase 0 corresponds to a noisy version of the true phase (p. Specifically, we consider the additive 
noise model 0 = (^ + e, where e is additive white Gaussian noise (AWGN). Given the wrapped image 
<0 = = W((^ + e), we would like to determine how accurately one can recover 0 in the presence of high 

levels of noise. Noisy phase images are typically more difficult to unwrap due to additional discontinuities 
that appear across the whole image. Note that robustness to noise is different from denoising, which implies 
that we do not attempt to reduce the level of noise during the unwrapping process. When unwrapping is 
successful, one can then denoise the unwrapped image with any state-of-the-art image denoising algorithm 
suitable to the noise at hand. 

In Table we report the SNR of the the estimate (f) with respect to the unwrapped phase (f) for GA, 
LS, PUMA, and IRTV algorithms. The true phase cp is the Man image of amplitude a = 6 and dimension 
128 X 128. The phase image <p is obtained by adding AWGN of variance corresponding to 16, 18, and 20 dB 
of input SNR (dB) = 10 log~ ^11^2)* table presents the median SNR for 10 independent 
realizations of the noise. In this example, the average running time of our MATLAB implementation of IRTV 
was about 7.5 minutes on our iMac using a 4 GHz Intel Core i7 processor. Additionally, the table presents 
the results of IRTV with different values of regularization parameter r. Specifically, we report the results 
for r G [0,10“^, 10“^, 10“^]. Higher levels of r imply stronger regularization during the reconstruction. The 
results in the table illustrate the advantage of using the Schatten norm of the Hessian to complement our 
data-fidelity term, i.e., one obtains a significant boost in unwrapping performance when r > 0. Moreover, 
the infiuence of r grows as the level of noise increases from 20 to 16 dB of input SNR. One must however 
note that similar to other regularization schemes, there is no theoretically optimal way of setting r and its 
optimal value might depend on the image, amount of wrapping, and noise. Our simulations indicate that 
the optimal value of r lies in the range [10“^, 10“^] for the configurations considered. 


4.2 Real Data 


In this second experimental part, we consider phase images that are acquired practically from distinct 
physical objects. In particular, we consider the setup for tomographic phase microscopy [GFYB+Q7 , which 
is a promising quantitative phase imaging technique. It is based on the principle that, for near-plane wave 
illumination of a sample, the phase of the transmitted field can be well approximated as the integral of the 
refractive index along the path of beam propagation. However, for the approximation to hold, the phases 
extracted from the transmitted fields must be first unwrapped, which significantly limits the applications of 
the technique for imaging objects that are off-focus, large, or have high index contrast. Once unwrapped, 
the phase image can simply be interpreted as the projection of refractive index, analogous to the projection 
of absorption in X-ray tomography. 

To evaluate our phase unwrapping algorithm, we measured refractive index tomograms of 6 polystyrene 
spheres (catalog no. 17135, Polysciences, refractive index n = 1.602 at 561 nm) immersed in oil with a lower 
refractive index of 1.516. The wrapped phase image of size 512 x 512 is extracted from the transmitted field 
at angle 39.02° with respect to vertical axis. This phase data is difficult to unwrap due of numerous visible 
phase discontinuities that appear along the borders of the beads. 

The parameters of our IRTV method were chosen as in the synthetic experiments above with the exception 
of the regularization parameter that was set to r = 10“^. Also as above, the reconstructed phase images 


were made consistent by using the operator (15). 

The results in Fig. [^illustrate the effectiveness of our method in unwrapping the phase even in the most 
difficult regions of the image that contain strong phase discontinuities. Specificaly, our method is the only 
one that was able to accurately unwrap the region between the two beads at the bottom of the image (see 
highlights in the figure). 
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5 Conclusion 


We have devised an algorithm for two-dimensional phase unwrapping that is based on an isotropic problem 
formulation. Based on suitable regularity assumptions, our technique has allowed to unwrap various phase 
images satisfactorily, including in the case where the phase contained significant amount of discontinuities. 
Compared to the standard techniques, the proposed method preserves edges of arbitrary orientations in the 
solution and effectively mitigates noise in practical configurations. From a general perspective, the obtained 
results further illustrate the interest of inverse-problem approaches for phase unwrapping. 
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